Ab initio GW many-body effects in graphene 
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We present an ab initio numerical many-body GW calculation of the band plot in free-standing 
graphene. We consider the full ionic and electronic structure introducing e-e interaction and cor- 
relation effects via a self-energy containing non-hermitian and dynamical terms. With respect to 
the density-functional theory local-density approximation, the Fermi velocity is renormalized with 
an increase of 17%, in better agreement with the experiment. Close to the Dirac point the linear 
dispersion is modified by the presence of a kink, as observed in angle-resolved photoemission spec- 
troscopy. We demonstrate that the kink is due to low-energy n — > n* single-particle excitations and 
to the 7r plasmon. The GW self-energy does not open the band gap. 

PACS numbers: 71.15.-m, 71.45.Gm, 79.20.Uv, 71.10.-w 
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The discovery of graphene by micromechanical cleav- 
age 0,0 and epitaxial grow [3j has attracted tremendous 
interest in consideration of its unusual electronic proper- 
ties. In the tight-binding (TB) formalism, the graphene 
2D honeycomb lattice structure gives rise to a semicon- 
ductor with zero band gap occurring at the K point in the 
Brillouin zone and a cone-like linear band-dispersion at 
low energy. This part is usually described by a massless 
Dirac (Weyl) dispersion 0]. Ab initio density- functional 
theory (DFT) calculations [j| confirm the TB linear dis- 
persion picture and give an estimate of the Fermi ve- 
locity vf lower by 15~20% than the experimental value. 
Recently, two angle-resolved photoemission spectroscopy 
(ARPES) experiments on graphene epitaxially grown on 
SiC S,J3, @, @, E3] raised the general interest. The first 
one [6j, la] observed at low energy a nearly linear band dis- 
persion with slight deviations in the form of small kinks 
interpreted as due to many-body electron-electron (e-e) 
and electronj)honon (e-ph) self-energy effects. The sec- 
ond one 0, H provided a different picture, with the 
opening of a band gap occurring at the Dirac K point and 
attributed either to substrate (SiC) or to many-body self- 
energy effects. A DFT calculation [llj seemed to confirm 
a substrate induced symmetry breaking, but recent STM 
measures 12j provided some evidence to exclude it. This 
situation calls for clarification about the role of e-e self- 
energy effects on the quasiparticle (QP) band plot, the 
Fermi velocity and the band gap opening. Previous ab 
initio works have dealt with e-ph effects [5|. |13I| and with 
e-e GW effects in graphene nanoribbons 1141 1. There are 
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also several non ab initio works 
studied e-e self-energy effects in a 2D massless Dirac 
model. 

In this work we calculate the band plot of free-standing 
undoped graphene introducing e-e interaction and cor- 
relation effects by an ab initio many-body GW self- 
energy We numerically simulate the full ionic 



and electronic structure of real graphene. We take 
into account the full dynamical dependence and non- 
hermiticity of the self-energy by an accurate contour- 
deformation (CD) integration. From the self-energy we 
then obtain the QP energies and the spectral function 
which can be directly compared with ARPES spectra. 
We show that the GW self-energy renormalizes the Fermi 
velocity by 17% such that it corrects the DFT underes- 
timation and leads to a value of 1.12 • 10 6 ms -1 , in good 
agreement with the accurate magnetotransport measure 
of 1.1 • 10 6 ms -1 2]. Furthermore, the nearly linear DFT 
band dispersion is in GW considerably distorted. Close 
to the Dirac point the self-energy results in an unusual 
negative GW band gap correction and the appearance of 
a kink in the band plot, leading to a scenario similar to 
that observed in ARPES, Ref. [6]. A comparison with the 
results of a GW plasmon-pole model calculation indicates 
that the kink is due to a coupling with the it plasmon at 
~ 5 eV and the low-energy it — > ir* single-particle excita- 
tions (SPE) shoulder present in the energy-loss function. 
This provides a partial confirmation to the explanation 
given in Ref. [f|. Finally, our results show that in free- 
standing graphene the GW self-energy does not open the 
band gap, in contrast to what is found in GW calcula- 
tions on graphene nanoribbons 14 1 . 



Our starting point is a standard ground-state DFT 
LDA calculation of infinite free-standing graphene. We 
use a plane waves basis set (62 Ry cutoff) and peri- 
odic boundary conditions on a hexagonal cell contain- 
ing 2 carbon atoms (a = 2.4 A) and 38 Bohr of vac- 
uum along the z direction, large enough to isolate spuri- 
ous replica of graphene sheets. We use Martins- Trouiller 
norm-conserving pseudopotentials with s and p electrons 
in the valence. We first calculate the ground-state energy 
and electronic density, and then the Kohn-Sham (KS) 
electronic structure to be used as starting point in the 
following GW excited-state calculation. 



2 



10- 



> 



CD 



-10 



-20- 









































** 



K 



M 



FIG. 1: Band plot of graphene. Solid thick lines: DFT-LDA 
KS; circles and dashed lines: GW. 



In the GW approximation ljj 2(| the self-energy is 
i I" 00 

£ GW (r, /, w) = — / dc^' Gfr, r', w - £«/)W(r, r', w') 

that is the product of the Green's function G and the dy- 
namically screened interaction W(uj) = £ _1 (w)w defined 
as the bare Coulombian interaction v screened by the 
dynamical dielectric function Vertex corrections 

are neglected both in the self-energy and in the polar- 
izability (hence in W). In the standard ab initio GW 
resolution procedure one builds G and W using the DFT 
KS electronic structure. The integral over the frequency 
is performed by a CD method [2l[ which consists in a de- 
formation of the real axis contour such as to calculate the 
self-energy as an integral along the imaginary axis minus 
a contribution arising from the residual of the contour- 
included poles of G. This is the most accurate method to 
perform a GW calculation. We also considered the stan- 
dard plasmon-pole model (PPM) approach [jo| . Once 
the integration is performed, we calculate the GW quasi- 
particle energies using a first-order perturbation theory 
expansion of £ around the LDA exchange-correlation po- 
tential v^ A and the KS energies ui — e KS 
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where Z = (1 — 9S /du>\ _ KS ) is the renormaliza- 
tion factor. We compare these energies to the positions 
of QP peaks in the spectral function 



A n (k,Lu) 
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calculated at the given k point in the Brillouin zone and 
projected on the considered n band. 

We used the ABINIT code. Convergence was achieved 
with 715 plane waves (10 Ry) to represent the wavefunc- 
tions and the exchange part T, x of the self-energy, 150 
and 200 bands to calculate respectively W and £ 27 1. 
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FIG. 2: Focus on the band plot Dirac point region. Solid thick 
line: DFT-LDA KS; circles and thin lines: GW CD; squares: 
GW PPM1 (u p ~ 15 eV); dashed line: fit over the GW CD 
linear region. 



The dimension of W and of the correlation part £ c as re- 
ciprocal space matrices was 169 (6 Ry). Imaginary axis 
integrations were carried out by a Gauss quadrature us- 
ing 10 nodes for the most accurate calculations and W 
was sampled every 0.2 eV over 100 frequencies on the real 
axis. The Brillouin zone was sampled with a (10 10 1) 
Monkhorst-Pack k-point grid. 

In Fig. Q] we compare the KS DFT-LDA (thick lines) 
and the quasiparticle GW (circles and dashed lines) 
electronic structures. Exchange end correlations effects 
slightly affect the band shapes. Relevant effects are a 
lowering of the a bands and an increase (up to +20%) of 
the gaps at M (4 -> 4.8 eV) and at T (6.4 -> 7 eV). This 
is a normal behavior of GW and in agreement with other 
calculations [H, [23| . We now focus on the Fermi energy 
Dirac point (K) region (Fig. [2]) . As previously obtained 
[|| , the DFT KS ir and n* (thick lines) band dispersion is 
linear in the first ~ 0.5 A -1 from the Dirac K point. The 
DFT KS Fermi velocity is 0.95 ■ 10 6 ms" 1 . This underes- 
timates by a 15% the experimental value. The dots and 
thin lines represent the GW band plot calculated by the 
CD method. The first evident self-energy effect is the loss 
of linearity along the region 0.05 A -1 close to the Dirac K 
point. However outside that region the linearity is soon 
recovered but with a slope larger than the DFT KS. A fit 
of the GW band with a straight line (dashed line) gives a 
Fermi velocity of 1 . 1 2 • 1 6 ms" 1 (1.14 with a th = 2.45 A). 
Thus the GW self-energy renormalizes by a +17% the 
DFT Fermi velocity and achieves a good agreement with 
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CD 


-0.023 


0.032 


-0.008 


0.729 +0.001 


0.328 


CD 


-0.011 


0.493 


-0.689 


~0.71 -0.121 


0.328 


PPM1 


-0.011 


0.493 


-0.316 


~ 0.76 +0.127 


0.328 


PPM2 


-0.011 


0.493 


-0.453 


0.537 +0.016 


0.328 


CD IRcut 


-0.011 


0.493 


-0.273 


0.771 +0.162 



TABLE I: Exchange and correlation components differences 
between bands 5 and 4, renormalization factor Z and GW 
band gap opening (energies in eV; Z as pure numbers). 



the experimental measure of 1.1 • 10 6 ms" 1 (2j. The resid- 
ual overestimation of 2~4% should be compensated by 
negative e-ph renormalization effects (—4% in Ref. fl3j|). 

We now focus on the non-linear region. In all the k- 
points far from the Dirac point, the GW correction acts 
in the usual direction to open the gap between DFT KS 
bands. On the other hand at k = (0.328,0.328,0) (re- 
duced coordinates), that is at ~ 0.025 A -1 from K, we 
have found an unusual negative GW correction of —0.12 
eV [28| which generates a kink at ~ 0.1 eV from the Dirac 
point. This result reproduces the experimental ARPES 
scenario of Ref. Q where a kink interpreted as due to e-e 
many-body effects is found more or less in this position. 
The position of our GW kink is also close to the position 
indicated in the other ARPES experiment (gray arrows 
in Fig. 3(d) of Ref. at 0.035 A" 1 and 0.17 eV). 

The negative GW correction conjuring the kink results 
from an unusual balance between the exchange and the 
correlation energy (see Tabled . Indeed, the exchange en- 
ergy difference between the bottom-of-conduction band 5 
and the top-of- valence band 4, (AY, X ), is typically several 
times larger than the difference in the correlation en ergy 
(AE C ) which opposes to the exchange energy [l9|, |2Cf . 
At the kink AE C is ~ 0.7 eV, larger than AS X ~ 0.5 eV. 
Therefore the negative correction to the band gap is due 
to a correlation energy stronger than usual. We report 
on Fig. [Hand Table Q] the result obtained by a PPM GW 
calculation (indicated as PPM1 and squares). Far from 
the Dirac point and exactly at the Dirac point (where the 
gap correction vanishes), the PPM GW bands precisely 
recover the CD result. PPM and CD start to deviate 
in the kink region and at the kink point the PPM pro- 
vides, in contrast to CD, a positive GW correction. This 
has an immediate interpretation. In the PPM approach 
the dynamical dependence of the energy-loss function 
— Se _1 entering into W and the GW self-energy, is repre- 
sented by a single plasmon-polc feature which is fitted to 
= 1 + Q 2 /(u 2 — a; 2 ). This is a good approximation 
in all systems where the energy-loss function presents a 
single plasmon feature. In graphene the energy-loss func- 
tion — 3e _1 (q ~ 0,w), as calculated from first-principles 
by the DP code in the RPA approximation (solid line 
in Fig. shows two major features: the total n + a 
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FIG. 3: Graphene energy-loss function at q ~ 0. Solid line: 
RPA with local- field effects (LFE); arrows: positions of the 
PPM1 (dot-dashed) and PPM2 (dashed) poles. 



plasmon at ~ 15 eV and the n plasmon at 
Fig. 1(e) of Ref. 
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5 eV (see 
Furthermore, at 



and Ref. [25j) 
the lowest energies there is a shoulder due to the it — > it* 
SPE. In our PPM1 GW calculation, the plasmon-pole fre- 
quency at q ~ sets to ui p = 15 eV (dot-dashed arrow) 
at the left side of the total plasmon peak. It reproduces 
the main energy-loss feature and implicitly accounts also 
for the low energy part of the spectrum. Other PPM 
calculations are examined (PPM2 in Fig. [3] and Table UJ). 
Forcing the plasmon-pole to adjust to lu p ~ 5 eV, the low 
energy part of the energy-loss is explicitly considered. In 
this case the resulting GW correction is around 0, closer 
to the correct GW CD result. Finally, we performed a 
CD calculation where e" 1 is computed cutting off low 
energy n — > n* SPE (2.5 eV IR cutoff). The low energy 
shoulder is suppressed (dotted line in Fig. [3]) and the 
intensity of the n plasmon is also unavoidably reduced. 
The IR cutoff has the effect to produce a positive GW 
correction even beyond the PPM1 result, thus validating 
the good quality of the PPM. From all these results we 
can deduce that both the tt plasmon and the low energy 
7r — > 7r* SPE shoulder provide the crucial contribution to 
the correlation energy and play a major role in conjuring 
the negative GW correction and the kink. 

The last important result of our GW calculation is that 
many-body effects, within the numerical error bar, do not 
open the band gap at the Dirac point. Both the CD and 
the PPM do not change the DFT-LDA band gap and 
band 4 and 5 keep degenerate at K. [2^] We find no indi- 
cation for a band gap opening even when looking at the 
spectral function. In Fig. 0] we reported the real (bottom 
panel) and the module of the imaginary part (middle) of 
the self-energy at the Dirac K point projected onto bands 
4 (solid line) and 5 (dashed line) . The intersection of the 
real part of H n k(u) with the straight line u> — e^J 5 + 
gives the position of the QP peak. Although KE is differ- 
ent for bands 4 and 5, the intersections occur at the same 
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FIG. 4: Graphene real and imaginary part of the self-energy 
and spectral function projected on bands 4 and 5 and at sev- 
eral k-points. 



point, so that the QP peaks are degenerate. The corre- 
spondent is in practice numerically 0. In the spectral 
function (top panel) we observe well defined QP peaks 
and weak 7r and ir + a plasmon satellites pointing to a 
normal Fermi liquid, unlike Rcf.s 15, HO, HI]. Thus 
our calculation excludes a band gap opening, even only 
apparent [181 ] , induced by a significant transfer of weigth 
from QP peaks to strong satellite plasmarons [lq . 

Finally we want to discuss possible mechanisms which 
can explain the presented picture. A simple kink, i.e. 
a slope change of the bands in one point, has been ob- 
tained in Ref. [26| by invoking a purely dynamical, k- 
independent, mechanism; this model is most adequate 
for strongly correlated systems like some transition metal 
oxides. The situation in graphene, with its negative GW 
correction and consequent s-shaped kink, requires a more 
complex description. We can indeed understand the re- 
sults by using a simplified PPM which accounts only for 
the low-lyin g tt — » 7r* SPE and the n plasmon. It has 
been shown [251 ] that for this plasmon we can assume a 
linear dispersion Lu p (q) = oj® + aq with a finite oj® ^ 0. 
Together with other straightforward assumptions, this 
model leads to a correlation energy in the form 



VfU>. 



q [log((«f + a)k + lu° p ) - log(afc + w°)] , 



which in the limit k — > oo tends to a constant whereas in 
the limit k — > fcj? is linear like the exchange energy ~E X ~ 
k, but with a coefficient VL 2 /2oj® times the exchange. 
This hypothesis is compatible with a zero band gap at 
k = kp and yields a large- fc region where the exchange 
dominates and renormalizes the Fermi velocity. If £1 > 
2u)p, it yields a region where the correlation dominates 
and provides negative GW corrections (the kink). 



In conclusion, we presented an ab initio many-body 
GW calculation in graphene. The GW self-energy renor- 
malizes by a 17% the DFT Fermi velocity but does not 
open a gap at K. Close to the Dirac K point the linear 
band dispersion is considerably affected by correlation, 
leading to the appearance of a kink. 
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Hartree-Fock and Coulomb-hole plus screened-exchange 
approximation [ll| calculations open the gap by 9 and 
5 meV respectively (34 and 13 meV at self-consistence, 
already achieved at the 2nd iteration) , all below the typ- 
ical ~ ±0.1 eV numerical error of ab initio calculations 



and so compatible with band gap. 



